Direct mapping of the finite temperature phase diagram of strongly correlated 

quantum models 
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Optical lattice experiments, with the unique potential of tuning interactions and density, have 
emerged as emulators of nontrivial theoretical models that are directly relevant for strongly corre- 
lated materials P, 0, 1, 1, 1, 1, EMM- 

However, so far the finite temperature phase diagram has not 
been mapped out for any strongly correlated quantum model. We propose a remarkable method for 
obtaining such a phase diagram for the first time directly from experiments using only the density 
profile in the trap as the input. We illustrate the procedure explicitly for the Bose Hubbard model, 
a textbook example of a quantum phase transition from a superfluid to a Mott insulator [10j]. Using 
"exact quantum Monte Carlo simulations in a trap with up to 10 6 bosons, we show that kinks in the 
local compressibility, arising from critical fluctuations, demarcate the boundaries between superfluid 
and normal phases in the trap. The temperature of the bosons in the optical lattice is determined 
from the density profile at the edge. Our method can be applied to other phase transitions even 
when reliable numerical results are not available. 



The grand challenge of condensed matter physics is to 
understand the emergence of novel phases arising from 
the organization of many degrees of freedom especially 
in regimes where particle interactions dominate over the 
kinetic energy. Brought to the forefront by the discov- 
ery of high temperature superconductivity in complex 
copper-based oxides, it is absolutely astounding that a 
simple model, like the Hubbard model is able to capture 
so many essential aspects of the physics of these mate- 
rials. However, whether the repulsive fermion Hubbard 
model really contains a d-wave superconducting ground 
state is still an open question after decades of study. 

Optical lattice experiments with cold atoms are emerg- 
ing as an amazing laboratory for making realizations of 
such bose and fermi Hubbard and Hciscnberg-type mod- 
els and observing phase transitions without the uncer- 
tainty posed by the complex materials. In this back- 
drop quantum Monte Carlo simulations in strongly in- 
teracting regimes are emerging as an important bridge 
between materials based condensed matter physics and 
cold atoms, highlighted in Fig. 1. The scale of the nu- 
merical simulations possible today with up to 10 6 parti- 
cles is able to match the experimental cold atom systems 
thereby allowing a direct comparison. 

The most important task for the quantum emulator 
of a Hamiltonian H is to determine the corresponding 
phase diagram in the temperature T , chemical potential 
fi and g space, where the latter are parameters of H. So 
far such a mapping has not been determined experimen- 
tally in any model. Even for the simplest Bose Hubbard 
model there have been many challenges arising from (i) 
lack of a clear diagnostic of how to identify phases. Tra- 
ditional methods using the sharpness of the interference 
patterns are not reliable to distinguish the normal and 
superfluid phases; (ii) complications due to coexistence 
of different phases in the same confining potential; (iii) 



FIG. 1: Interaction between novel materials, models and op- 
tical lattice emulators. The models capture the fundamental 
physics of complex materials with novel phenomena such as 
high temperature superconductivity and quantum magnetism 
and form a bridge with experiments on ultracold atoms in 
optical lattices that emulate these models in a clean environ- 
ment. 



lack of thermometry of the Bose gas in the optical lattice; 
(iv) and until very recently, the lack of QMC simulations 
on experimentally relevant sizes. The work reported here 
overcomes all of the above obstacles. 

The Hamiltonian for the single band Bose Hubbard 
Model (BHM) is given by: 
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where bi(b\) is the boson destruction (creation) oper- 
ator at a site i, rii — b\bi is the density operator, z = 6 
is the coordination number in 3D, and /io controls the 
density of bosons. The relative strength, g = t/U of 
the tunnelling t of bosons between nearest neighbor sites 
vs. the repulsive interaction U between bosons, tunes 
the system through a superfluid to Mott transition at 
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FIG. 2: Phase diagram of BHM. (a) fi — t/U plane at T = (schematic); (b) fi — t/U plane at T (schematic). The superfluid 
phase is characterised by a non-zero order parameter (6^) ^ and a non-zero superfluid density p s 7^ 0, both of which vanish 
in the normal phase. The Mott insulator at T = 0, characterized by integer filling and a finite gap to excitations Ao, crosses 
over to the normal state for T ~ Ao . At finite temperatures, the critical chemical potential p, c (T,t/U) demarcates the S-N 
phase boundary. In trapped atomic gases, the local chemical potential decreases from the centre of the trap to the edge (brown 
vertical line), (c) QMC simulations and finite size scaling to calculate the phase diagram in the p. — T plane for a fixed t/U. 



T = (Fig. 2(a)). At finite T the system shows a phase 
transition from a superfluid (S) to a normal(N) phase 
(Fig. 2(b)). In optical lattice experiments, bosons are 
confined in an additional harmonic trap of frequency ui 



modelled by H = Hbh m 



where m is the 



mass of the bosons. We simulate the Hamiltonian H with 
the trap using quantum Monte Carlo techniques that in- 
clude the effect of strong interactions "exactly" within 
statistical errors. Recent modifications of the directed- 
loop al gori thm for the world-line quantum Monte Carlo 
method [111] have allowed us to significantly improve the 
efficiency near a critical point [TJ] for large 3D systems 
with up to 10 6 bosons in a 64 3 lattice for the first time. 
Thus the phase diagram in the p. — T plane is calculated 
by finite size scaling for fixed t/U .(Fig. 2(c)). 

Experimentally, the challenge of obtaining the phase 
diagram at finite T lies in identifying measurable prop- 
erties that can diagnose different phases. Tradition- 
ally the momentum distribution n(k) imaged in the in- 
terference patterns of the expanding cold atom clouds 
has been used to identify the phases. In a long- 
time ballistic expansion the interference pattern fi{f) = 

(m/hrf \W (k = ) \ 2 n (k = ff ) essentially provides 

an image of n(k) — ^2 ii j(b i bj)e lk '^ n ~ r: >^ before expan- 
sion by convoluting with W(k) , the Fourier transform 
of the Wannier function within a single site. Here r is 
the expansion time. The final image detects the col- 
umn integrated momentum distribution Nj_(k x ,k y ) = 
J dk z \W(k)\ 2 n(k). We have previously shown for a ho- 
mogeneous BHM that sharp peaks in the interference 
pattern are not reliable for identifying a superfluid 13 . 
Ibl . Il5j |. What is the effect of a confining potential on 
n(jfe)? 

To answer this question we calculate both n(k) and 
the density profile p(r) in a harmonic trap at different T 
and tuning t/U (Fig. 3). For comparison, within local 



density approximation (LDA), we also calculate p h (p(r)), 
the density p h (p) for a homogeneous system using QMC, 
where p,(r) = p,o — muj 2 r 2 /I. Fig. 3 shows the excellent 
agreement of p(r) and p h (p(r)), not unexpected for a 
shallow harmonic trap. We next turn to a key diagnos- 
tic, the local superfluid density distribution in the trap 
Ps (/ i ( r ))i obtained from a knowledge of fi(r) and the su- 
perfluid density Pg(p>) in the homogeneous system. The 
most significant observation from these data is that even 
in the presence of a trap when all the atoms are in the 
normal state, the interference pattern continues to show 
sharp peaks. Strong repulsive interactions that expand 
the bosonic cloud and suppress T c in the optical lattice 
contribute to a sharpening up of n(k) even in the normal 
state. With the lowering of T, as more regions within the 
trap become superfluid, there is a distinct change in the 
shape of the sharp peak (Fig. 3(d)). The emergence of 
a singular feature of width ~ 1/L, limited by the cloud 
size, indicates phase coherence throughout the sample. 

It is evident from Fig. 3 that n(k) which integrates 
over the entire trap is unable to provide local informa- 
tion about the distribution of the phases. Also the finite 
expansion time, as well as resolution problems compli- 
cates the analysis of the interference pattern h(r) and its 
relation to n(fc)[l6j]. We therefore explore direct methods 
of identifying phase boundaries in the trap. 

Our proposal relies on extracting the local compress- 
ibility in the trap defined by Kdiff 



from a 

high- resolution scan of p(r) in a trap. Kdiff agrees very 
well with the compressibility K h (p(r)) = for a bulk 
system (Fig. 4(a)). A novel feature of the local compress- 
ibility is the existence of kinks at specific locations in 
the trap. The origin for the kink feature becomes evi- 
dent if we consider in addition the behavior of the local 
superfluid density in the trap p^(p(r)) (Fig. 5). We see 
clearly that the kinks in the compressibility occur at the 
S-N phase boundary. At fixed T, the superfluid order 
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FIG. 3: Interference pattern in a trap with TV w 10 bosons: 
The insets show the local density profile p(r)(red boxes) ob- 
tained using QMC as a function of the radial coordinate in 
the trap. Also shown is the excellent agreement with p h (p(r)) 
calculated within local density approximation (LDA) (purple 
diamonds). The superfluid density p^(p(r)) is non-zero in a 
significant portion of the trap only at sufficiently low temper- 
atures in panel (d). Note that sharp peaks are seen in n(k) 
in panels (a), (b) and (c) even when all the atoms are in the 
normal state. The peak in panel (d) is much sharper. The 
difference between the normal state and the superfluid state 
is evident from the system-size dependence of the width of 
the peak. 



becomes weaker when approaching the phase boundary. 
The increase of fluctuations is reflected in the growth of 
K h (r). Starting from deep in the Mott-like region, where 
the particle number per site is essentially unity, n h (r) 
increases when approaching the phase boundary with a 
superfluid resulting in a singular feature in K h (r). The 
phase transition in the BHM is in the XY universality 
class with n h (r) ~ |/x— p. c \~ a and a = — 1 giving a cusp at 
the transition. The singular feature is evident even with 
rounding from finite size effects in a trap as seen from 
Fig 4 and 5. At high temperatures, when p s vanishes in 
the trap, the kinks disappear in the absence of a phase 
boundary and the compressibility changes smoothly. 

The importance of Fig. 4 and 5 is that purely from 
a high resolution scan of the local density it is possible 
to obtain the local compressibility and from the exis- 
tence of kinks deduce that there must be superfluid re- 
gions in the trap separated from normal regions, even 
though both of them have non-integer filling and the den- 
sity changes smoothly across the phase boundary. We 
would like to point out that there are significant differ- 
ences between our work and previous studies that have 
primarily focused on the compressibility of Mott and su- 



18l . [19|. As seen from Fig. 



pcrfluid phases at T 
5(a), the compressibility of the Mott state is exactly zero 
and becomes finite in the superfluid. At finite T, noth- 
ing singular happens as we move from Mott-like to more 
compressible regions in the trap. It is only at the S- 
N boundary (Fig. 5(b,c,d)) even though formed between 
two compressible phases, that the system shows enhanced 
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FIG. 4: Kinks in the local compressibility and mapping the 
phase diagram: (a) The local density p(r)(blue dots) agrees 
with p h (r) (purple diamonds) calculated for a homogeneous 
system. We also show that the compressibility Kdiff{i") (red 
boxes) directly obtained from the density profile agrees with 
K h (p(r)) (pink triangles) calculated within LDA. The signal 
to noise in the derivative method to extract Kdiff(r) can be 
improved by taking an angular average of p(r). Experimen- 
tally, p(r) is obtained by an inverse Abel transformation of 
the column density[20J or by identifying the planar-integrated 
density with the pressure. Recently, the column density has 
been measured in high resolution experiments [2l[. Alterna- 
tively, p(x,y, Zq) at fixed zq can be obtained directly [2^| using 
a sliding technique to get the density, (b) Determination of 
temperature by fitting (shown by lines) of the tail of the den- 
sity profile (QMC results shown by symbols) to the ideal gas 
behavior in the dilute regime e 13 ^ <C 1 and \p(r) + 2t\ 3> Un so 
interactions can be ignored. Inset shows the fitting to the an- 
gular integrated densities at the corresponding temperatures 
that have larger signals and therefore better accuracy. T/t 
are 0.67, 0.57 and 0.44 from top to bottom. 



critical fluctuations detected in the local compressibility. 
Note that in Fig. 5(b, d) there are kinks even in the 
absence of any Mott regions. 

Usually strong interaction effects in an optical lat- 
tice hamper the determination of the temperature T 
and the chemical potential (1q at the center in the op- 
tical lattice. We propose two methods: (i) By fitting 
p(r) with p(/i(r),T) constrained by the known number 
of bosons in the trap according to N = Jd 3 rp(r) = 

Ho ^ fy ^). (ii) By fitting the tail of the den- 
sity profile in the trap to the expected ideal gas behav- 

iOT Pldeal = J (>(e fe -ft,+m"V/2) _ {j- 1 as demon _ 

strated in Fig. 4(b). This is a direct method with- 
out any input from simulations. Thus from the location 
r c of kinks in Kdiff{r), experimentalists can determine 
/i c = Po — mu! 2 r^/2 to map out the phase diagram in the 
fJ> — T plane for a fixed t/U and compare with theoretical 
results in Fig. 2(c). 

Mapping the phase diagrams of strongly correlated 
quantum models is the central goal of condensed matter 
physics. In that context our proposed method of prob- 
ing fluctuations of the density, specifically singularites or 
cusps, arising from critical fluctuations, provides a gen- 
eral method for identifying the phase boundaries. The 
cold atom emulator of a mathematical model in a confin- 
ing trap, essentially generates a chemical potential scan 
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FIG. 5: Origin of kinks: The local density p(/i(r)) (purple), 
compressibility re h (/i(r))(red) and superfluid density 
ps (/x(r))(blue) as functions of the radial coordinate in 
the trap. Panel (a) shows when p a — in the trap, /c(r)shows 
a smooth variation. As T is reduced or t/U is increased 
(b,c,d) and a finite p s develops in some portion of the trap, 
fc(r) shows sharp kinks. The location of these kinks coincides 
with the S-N boundary in the trap. Notice the location of 
the kinks shifts as the S-N boundary changes. 



of the phase diagram from a single measurement. Our 
studies open up new directions for mapping the finite 
temperature phase diagrams of strongly correlated quan- 
tum models, usually not available by other means, by 
probing local properties of trapped atoms in optical lat- 
tices rather than relying on the momentum distribution. 
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